resultsPath=file.path(getwd(),"Results")
# Gather parameters from command line
#dir.create(file.path(resultsPath,"cache"), showWarnings=F, recursive=T)
nCores <- parallel::detectCores()#params$nCores
subsetGenes <- params$subsetGenes
subsetCells <- params$subsetCells
resolution <- as.numeric(params$resolution)
root <- getwd()
# Have to setwd via knitr
# knitr::opts_knit$set(root.dir=resultsPath, child.path = resultsPath)
knitr::opts_chunk$set(echo=T, error=T, root.dir = resultsPath#cache=T, cache.lazy=T,
)
# kableStyle = c("striped", "hover", "condensed", "responsive")
# Utilize parallel processing later on
print(paste("**** Utilized Cores **** =", nCores)) ## [1] "**** Utilized Cores **** = 4"
params## $resultsPath
## [1] "./"
##
## $subsetGenes
## [1] "protein_coding"
##
## $subsetCells
## [1] 500
##
## $resolution
## [1] 0.6
** ./ **
library(Seurat)
library(dplyr)
library(gridExtra)
library(knitr)
library(plotly)
library(ggplot2)
library(reshape2)
library(shiny)
library(ggrepel)
library(DT)
#
# install.packages('devtools')
# devtools::install_github('talgalili/heatmaply')
## Install Bioconductor
# if (!requireNamespace("BiocManager"))
# install.packages("BiocManager")
# BiocManager::install(c("biomaRt"))
library(biomaRt)
# BiocManager::install(c("DESeq2"))
library(DESeq2)
createDT <- function(DF, caption="", scrollY=500){
data <- datatable(DF, caption=caption,
extensions = 'Buttons',
options = list( dom = 'Bfrtip',
buttons = c('copy', 'csv', 'excel', 'pdf', 'print'),
scrollY = scrollY, paging = F
)
)
return(data)
}
# Useful Seurat functions
## Seurat::FindGeneTerms() # Enrichr API
## Seurat::MultiModal_CCA() # Integrates data from disparate datasets (CIA version too)
sessionInfo()## R version 3.5.1 (2018-07-02)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS 10.14.2
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] DESeq2_1.22.2 SummarizedExperiment_1.12.0
## [3] DelayedArray_0.8.0 BiocParallel_1.16.5
## [5] matrixStats_0.54.0 Biobase_2.42.0
## [7] GenomicRanges_1.34.0 GenomeInfoDb_1.18.1
## [9] IRanges_2.16.0 S4Vectors_0.20.1
## [11] BiocGenerics_0.28.0 biomaRt_2.38.0
## [13] DT_0.5.1 ggrepel_0.8.0
## [15] shiny_1.2.0 reshape2_1.4.3
## [17] plotly_4.8.0 knitr_1.21
## [19] gridExtra_2.3 dplyr_0.7.8
## [21] Seurat_2.3.4 Matrix_1.2-15
## [23] cowplot_0.9.3 ggplot2_3.1.0
##
## loaded via a namespace (and not attached):
## [1] snow_0.4-3 backports_1.1.3 Hmisc_4.1-1
## [4] plyr_1.8.4 igraph_1.2.2 lazyeval_0.2.1
## [7] splines_3.5.1 digest_0.6.18 foreach_1.4.4
## [10] htmltools_0.3.6 lars_1.2 gdata_2.18.0
## [13] magrittr_1.5 checkmate_1.8.5 memoise_1.1.0
## [16] cluster_2.0.7-1 mixtools_1.1.0 ROCR_1.0-7
## [19] annotate_1.60.0 R.utils_2.7.0 prettyunits_1.0.2
## [22] colorspace_1.3-2 blob_1.1.1 xfun_0.4
## [25] crayon_1.3.4 RCurl_1.95-4.11 jsonlite_1.6
## [28] genefilter_1.64.0 bindr_0.1.1 survival_2.43-3
## [31] zoo_1.8-4 iterators_1.0.10 ape_5.2
## [34] glue_1.3.0 gtable_0.2.0 zlibbioc_1.28.0
## [37] XVector_0.22.0 kernlab_0.9-27 prabclus_2.2-6
## [40] DEoptimR_1.0-8 scales_1.0.0 mvtnorm_1.0-8
## [43] DBI_1.0.0 bibtex_0.4.2 Rcpp_1.0.0
## [46] metap_1.0 dtw_1.20-1 viridisLite_0.3.0
## [49] xtable_1.8-3 progress_1.2.0 htmlTable_1.13.1
## [52] reticulate_1.10 foreign_0.8-71 bit_1.1-14
## [55] proxy_0.4-22 mclust_5.4.2 SDMTools_1.1-221
## [58] Formula_1.2-3 tsne_0.1-3 htmlwidgets_1.3
## [61] httr_1.4.0 gplots_3.0.1 RColorBrewer_1.1-2
## [64] fpc_2.1-11.1 acepack_1.4.1 modeltools_0.2-22
## [67] ica_1.0-2 pkgconfig_2.0.2 XML_3.98-1.16
## [70] R.methodsS3_1.7.1 flexmix_2.3-14 nnet_7.3-12
## [73] locfit_1.5-9.1 tidyselect_0.2.5 rlang_0.3.0.1
## [76] later_0.7.5 AnnotationDbi_1.44.0 munsell_0.5.0
## [79] tools_3.5.1 RSQLite_2.1.1 ggridges_0.5.1
## [82] evaluate_0.12 stringr_1.3.1 yaml_2.2.0
## [85] npsurv_0.4-0 bit64_0.9-7 fitdistrplus_1.0-11
## [88] robustbase_0.93-3 caTools_1.17.1.1 purrr_0.2.5
## [91] RANN_2.6 bindrcpp_0.2.2 pbapply_1.3-4
## [94] nlme_3.1-137 mime_0.6 R.oo_1.22.0
## [97] hdf5r_1.0.1 compiler_3.5.1 rstudioapi_0.8
## [100] png_0.1-7 lsei_1.2-0 geneplotter_1.60.0
## [103] tibble_2.0.0 stringi_1.2.4 lattice_0.20-38
## [106] trimcluster_0.1-2.1 pillar_1.3.1 Rdpack_0.10-1
## [109] lmtest_0.9-36 data.table_1.11.8 bitops_1.0-6
## [112] irlba_2.3.2 gbRd_0.4-11 httpuv_1.4.5.1
## [115] R6_2.3.0 latticeExtra_0.6-28 promises_1.0.1
## [118] KernSmooth_2.23-15 codetools_0.2-16 MASS_7.3-51.1
## [121] gtools_3.8.1 assertthat_0.2.0 withr_2.1.2
## [124] GenomeInfoDbData_1.2.0 diptest_0.75-7 doSNOW_1.0.16
## [127] hms_0.4.2 grid_3.5.1 rpart_4.1-13
## [130] tidyr_0.8.2 class_7.3-15 rmarkdown_1.11
## [133] segmented_0.5-3.0 Rtsne_0.15 base64enc_0.1-3
print(paste("Seurat ", packageVersion("Seurat")))## [1] "Seurat 2.3.4"
setwd("~/Desktop/PD_scRNAseq/")
dir.create(file.path(root,"Data"), showWarnings=F)
load(file.path(root,"Data/seurat_object_add_HTO_ids.Rdata"))
pbmc <- seurat.obj
rm(seurat.obj)pbmc## An object of class seurat in project RAJ_13357
## 24914 genes across 22113 samples.
metadata <- read.table(file.path(root,"Data/meta.data4.tsv"))
createDT( metadata, caption = "Metadata") ## Warning in instance$preRenderHook(instance): It seems your data is too
## big for client-side DataTables. You may consider server-side processing:
## https://rstudio.github.io/DT/server.html
# Make AgeGroups
makeAgeGroups <- function(){
dim(metadata)
getMaxRound <- function(vals=metadata$Age, unit=10)unit*ceiling((max(vals)/unit))
getMinRound <- function(vals=metadata$Age, unit=10)unit*floor((min(vals)/unit))
ageBreaks = c(seq(getMinRound(), getMaxRound(), by = 10), getMaxRound()+10)
AgeGroupsUniq <- c()
for (i in 1:(length(ageBreaks)-1)){
AgeGroupsUniq <- append(AgeGroupsUniq, paste(ageBreaks[i],ageBreaks[i+1], sep="-"))
}
data.table::setDT(metadata,keep.rownames = T,check.names = F)[, AgeGroups := cut(Age,
breaks = ageBreaks,
right = F,
labels = AgeGroupsUniq,
nclude.lowest=T)]
metadata <- data.frame(metadata)
unique(metadata$AgeGroups)
head(metadata)
dim(metadata)
return(metadata)
}
# metadata <- makeAgeGroups()
pbmc <- AddMetaData(object = pbmc, metadata = metadata)
# Get rid of any NAs (cells that don't match up with the metadata)
if(subsetCells==F){
pbmc <- FilterCells(object = pbmc, subset.names = "nGene", low.thresholds = 0)
} else {pbmc <- FilterCells(object = pbmc, subset.names = "nGene", low.thresholds = 0,
# Subset for testing
cells.use = pbmc@cell.names[0:subsetCells]
)
} Include only subsets of genes by type. Biotypes from: https://useast.ensembl.org/info/genome/genebuild/biotypes.html
subsetBiotypes <- function(pbmc, subsetGenes){
if( subsetGenes!=F ){
print(paste("Subsetting genes:",subsetGenes))
# If the gene_biotypes file exists, import csv. Otherwise, get from biomaRt
if(file_test("-f", file.path(root,"Data/gene_biotypes.csv"))){
biotypes <- read.csv(file.path(root,"Data/gene_biotypes.csv"))
}
else {
ensembl <- useMart(biomart="ENSEMBL_MART_ENSEMBL", host="grch37.ensembl.org",
dataset="hsapiens_gene_ensembl")
ensembl <- useDataset(mart = ensembl, dataset = "hsapiens_gene_ensembl")
listFilters(ensembl)
listAttributes(ensembl)
biotypes <- getBM(attributes=c("hgnc_symbol", "gene_biotype"), filters="hgnc_symbol",
values=row.names(pbmc@data), mart=ensembl)
write.csv(biotypes, file.path(root,"Data/gene_biotypes.csv"), quote=F, row.names=F)
}
# Subset data by creating new Seurat object (annoying but necessary)
geneSubset <- biotypes[biotypes$gene_biotype==subsetGenes,"hgnc_symbol"]
print(paste(dim(pbmc@raw.data[geneSubset, ])[1],"/", dim(pbmc@raw.data)[1],
"genes are", subsetGenes))
# Add back into pbmc
subset.matrix <- pbmc@raw.data[geneSubset, ] # Pull the raw expression matrix from the original Seurat object containing only the genes of interest
pbmc_sub <- CreateSeuratObject(subset.matrix) # Create a new Seurat object with just the genes of interest
orig.ident <- row.names(pbmc@meta.data) # Pull the identities from the original Seurat object as a data.frame
pbmc_sub <- AddMetaData(object = pbmc_sub, metadata = pbmc@meta.data) # Add the idents to the meta.data slot
pbmc_sub <- SetAllIdent(object = pbmc_sub, id = "ident") # Assign identities for the new Seurat object
pbmc <- pbmc_sub
rm(list = c("pbmc_sub","geneSubset", "subset.matrix", "orig.ident"))
}
}
subsetBiotypes(pbmc, subsetGenes)## [1] "Subsetting genes: protein_coding"
## [1] "14827 / 24914 genes are protein_coding"
Filter by cells, normalize , filter by gene variability.
pbmc <- FilterCells(object = pbmc, subset.names = c("nGene", "percent.mito"),
low.thresholds = c(200, -Inf), high.thresholds = c(2500, 0.05))
pbmc <- NormalizeData(object = pbmc, normalization.method = "LogNormalize",
scale.factor = 10000)** Important!**: Specify do.par = T, and num.cores = nCores in ‘ScaleData’ to use all available cores.
# Store the top most variable genes in @var.genes
pbmc <- FindVariableGenes(object = pbmc, mean.function = ExpMean, dispersion.function = LogVMR,
x.low.cutoff = 0.0125, x.high.cutoff = 3, y.cutoff = 0.5)# IMPORTANT!: Must set do.par=T and num.cors = n for large datasets being processed on computing clusters
pbmc <- ScaleData(object = pbmc, vars.to.regress = c("nUMI", "percent.mito"), do.par = T, num.cores = nCores)## Regressing out: nUMI, percent.mito
##
## Time Elapsed: 11.3517928123474 secs
## Scaling data matrix
pbmc## An object of class seurat in project RAJ_13357
## 24914 genes across 495 samples.
vp <- VlnPlot(object = pbmc, features.plot = c("nGene", "nUMI", "percent.mito"),nCol = 3, do.return = T) %>% + ggplot2::aes(alpha=0.5)
vp# par(mfrow = c(1, 2))
gp1 <- GenePlot(object = pbmc, gene1 = "nUMI", gene2 = "percent.mito", pch.use=20,
do.hover=T, data.hover = "mut")gp1gp2 <- GenePlot(object = pbmc, gene1 = "nUMI", gene2 = "nGene", pch.use=20,
do.hover=T, data.hover = "mut")gp2ProjectPCA scores each gene in the dataset (including genes not included in the PCA) based on their correlation with the calculated components. Though we don’t use this further here, it can be used to identify markers that are strongly correlated with cellular heterogeneity, but may not have passed through variable gene selection. The results of the projected PCA can be explored by setting use.full=T in the functions above
# Run PCA with only the top most variables genes
pbmc <- RunPCA(object = pbmc, pc.genes = pbmc@var.genes, do.print=F) #, pcs.print = 1:5, genes.print = 5VizPCA(object = pbmc, pcs.use = 1:2)PCAPlot(object = pbmc, dim.1 = 1, dim.2 = 2, do.hover=T, data.hover="mut")pbmc <- ProjectPCA(object = pbmc, do.print=F)
# 'PCHeatmap' is a wrapper for heatmap.2
# PCA Heatmap: PC1-PCn
PCHeatmap(object = pbmc, pc.use = 1:12, do.balanced=T, label.columns=F, use.full=F) #
# PCHeatmap_interactive <- function(PC=1){
# PC_dat <- PCHeatmap(object = pbmc, pc.use = PC, do.return = T)
# # Cluster samples
# Xclust <- pcp %>% dist(upper = T) %>% hclust()
# Yclust <- PC_dat %>% t() %>% dist(upper = T) %>% hclust()
# PC_dat <- PC_dat[Xclust$order, Yclust$order]
# # Plotly
# plot_ly(y=row.names(PC_dat), z=matrix(PC_dat), type = "heatmap",
# colors =viridis::plasma(n=100))
# }
# PCHeatmap_interactive(PC=1) Determine statistically significant PCs for further analysis. NOTE: This process can take a long time for big datasets, comment out for expediency. More approximate techniques such as those implemented in PCElbowPlot() can be used to reduce computation time
#pbmc <- JackStraw(object = pbmc, num.replicate = 100, display.progress = FALSE)
PCElbowPlot(object = pbmc)We first construct a KNN graph based on the euclidean distance in PCA space, and refine the edge weights between any two cells based on the shared overlap in their local neighborhoods (Jaccard similarity). To cluster the cells, we apply modularity optimization techniques such as the Louvain algorithm (default) or SLM [SLM, Blondel et al., Journal of Statistical Mechanics], to iteratively group cells together, with the goal of optimizing the standard modularity function.
On Resolution
The FindClusters function implements the procedure, and contains a resolution parameter that sets the ‘granularity’ of the downstream clustering, with increased values leading to a greater number of clusters. We find that setting this parameter between 0.6-1.2 typically returns good results for single cell datasets of around 3K cells. Optimal resolution often increases for larger datasets. The clusters are saved in the object@ident slot.
# TRY DIFFERENT RESOLUTIONS
pbmc <- StashIdent(object = pbmc, save.name = "pre_clustering")
# pbmc <- SetAllIdent(object = pbmc, id = "pre_clustering")
pbmc <- FindClusters(object = pbmc, reduction.type = "pca", dims.use = 1:10,
resolution = resolution, print.output = 0, save.SNN = T)
PrintFindClustersParams(object = pbmc) ## Parameters used in latest FindClusters calculation run on: 2019-01-08 17:05:57
## =============================================================================
## Resolution: 0.6
## -----------------------------------------------------------------------------
## Modularity Function Algorithm n.start n.iter
## 1 1 100 10
## -----------------------------------------------------------------------------
## Reduction used k.param prune.SNN
## pca 30 0.0667
## -----------------------------------------------------------------------------
## Dims used in calculation
## =============================================================================
## 1 2 3 4 5 6 7 8 9 10
pbmc <- StashIdent(object = pbmc, save.name = "post_clustering") pbmc <- RunUMAP(object = pbmc, dims.use = 1:10)
# Plot results
DimPlot(object = pbmc, reduction.use = 'umap')As input to the tSNE, we suggest using the same PCs as input to the clustering analysis, although computing the tSNE based on scaled gene expression is also supported using the genes.use argument.
** Important!**: Specify num_threads=0 in ‘RunTSNE’ to use all available cores.
labSize <- 6
pbmc <- RunTSNE(object=pbmc, reduction.use = "pca", dims.use = 1:10, do.fast = TRUE,
tsne.method = "Rtsne", num_threads=0) # num_threads
# note that you can set do.label=T to help label individual clusters
TSNEPlot(object = pbmc, do.label=T, label.size = labSize, do.return=T) %>% ggplotly() tSNE_metadata_plot <- function(var){
print(paste("t-SNE Metadata plot for ", var))
# Metadata plot
p1 <- TSNEPlot(pbmc, do.return = T, do.label = T, group.by = var, pt.size=1,
plot.title=paste("Color by ",var), vector.friendly=T) %>% ggplotly() %>%
layout(legend = list(orientation = 'h', xanchor = "center", x = 0.5, y = .999))
# t-SNE clusters plot
p2 <- TSNEPlot(pbmc, do.return = T, do.label = T, pt.size=1,
plot.title=paste("Color by Clusters"), vector.friendly=T) %>% ggplotly() %>%
layout(legend = list(orientation = 'h', xanchor = "center", x = 0.5, y = .999))
#print(plot_grid(ggplotly(p1), ggplotly(p2)))
fluidPage(
fluidRow(
column(6, p1), column(6, p2)
)
)
}
# metaVars <- c(dx","mut","Gender","Age")
#
# for (var in metaVars){
# print(paste("t-SNE Metadata plot for ",var))
# # Metadata plot
# p1 <- TSNEPlot(pbmc, do.return = T, pt.size = 0.5, group.by = var, do.label = T,
# dark.theme=F, plot.title=paste("Color by ",var))
# # t-SNE clusters plot
# p2 <- TSNEPlot(pbmc, do.label = T, do.return = T, pt.size = 0.5, plot.title=paste("Color by t-SNE clusters"))
# print(plot_grid(p1, p2))
# } tSNE_metadata_plot("dx") ## [1] "t-SNE Metadata plot for dx"
tSNE_metadata_plot("mut") ## [1] "t-SNE Metadata plot for mut"
tSNE_metadata_plot("Gender") ## [1] "t-SNE Metadata plot for Gender"
tSNE_metadata_plot("Age") ## [1] "t-SNE Metadata plot for Age"
Seurat has several tests for differential expression which can be set with the test.use parameter (see the DE vignette for details). For example, the ROC test returns the ‘classification power’ for any individual marker (ranging from 0 - random, to 1 - perfect).
Shown here: Biomarkers of each cluster vs. all other clusters.
pbmc.markers <- FindAllMarkers(object = pbmc, only.pos = TRUE,
min.pct = 0.25, thresh.use = 0.25)
createDT(pbmc.markers, caption = paste("All Biomarkers: All Clusters"))topNum = 5
topBiomarkers <- pbmc.markers %>% group_by(cluster) %>% top_n(topNum, avg_logFC)
createDT(pbmc.markers, caption = paste("All Biomarkers: All Clusters"))getTopBiomarker <- function(pbmc.markers, clusterID, topN=1){
df <-pbmc.markers %>%
subset(p_val_adj<0.05 & cluster==as.character(clusterID)) %>%
arrange(desc(avg_logFC))
top_pct_markers <- df[1:topN,"gene"]
return(top_pct_markers)
}
# clust1_biomarkers <- getTopBiomarker(pbmc.markers, clusterID=1, topN=2)
# clust2_biomarkers <- getTopBiomarker(pbmc.markers, clusterID=2, topN=2)
### Plot biomarkers
plotBiomarkers <- function(pbmc, biomarkers, cluster){
biomarkerPlots <- list()
for (marker in biomarkers){
p <- VlnPlot(object = pbmc, features.plot = c(marker), y.log=T, return.plotlist=T)
biomarkerPlots[[marker]] <- p + ggplot2::aes(alpha=0.5) + xlab( "Cluster") + ylab( "Expression")
}
combinedPlot <- do.call(grid.arrange, c(biomarkerPlots, list(ncol=2, top=paste("Top DEG Biomarkers for Cluster",cluster))) )
# biomarkerPlots <- lapply(biomarkers, function(marker) {
# VlnPlot(object = pbmc, features.plot = c(marker), y.log=T, return.plotlist=T) %>% + ggplot2::ggtitle(marker) %>% ggplotly()
# })
# return(subplot(biomarkerPlots) )
}
top1 <- pbmc.markers %>% group_by(cluster) %>% top_n(1, avg_logFC)
nCols <- floor( sqrt(length(unique(top1$cluster))) )
figHeight <- nCols *7
# Plot top 2 biomarker genes for each
for (clust in unique(pbmc.markers$cluster)){
cat('\n')
cat("### Cluster ",clust,"\n")
biomarkers <- getTopBiomarker(pbmc.markers, clusterID=clust, topN=2)
plotBiomarkers(pbmc, biomarkers, clust)
cat('\n')
} ##Construct the plot object
volcanoPlot <- function(DEG_df, caption="", topFC_labeled=5){
DEG_df$sig<- ifelse( DEG_df$p_val_adj<0.05 & DEG_df$avg_logFC<1.5, "p_val_adj<0.05",
ifelse( DEG_df$p_val_adj<0.05 & DEG_df$avg_logFC>1.5, "p_val_adj<0.05 & avg_logFC>1.5",
"p_val_adj>0.05"
))
DEG_df <- arrange(DEG_df, desc(sig))
vol <- ggplot(data=DEG_df, aes(x=avg_logFC, y= -log10(p_val_adj))) +
geom_point(alpha=0.5, size=3, aes(col=sig)) +
scale_color_manual(values=list("p_val_adj<0.05"="turquoise3",
"p_val_adj<0.05 & avg_logFC>1.5"="purple",
"p_val_adj>0.05" = "darkgray")) +
theme(legend.position = "none") +
xlab(expression(paste("Average ",log^{2},"(fold change)"))) +
ylab(expression(paste(-log^{10},"(p-value)"))) + xlim(-2,2) +
## ggrepl labels
geom_text_repel(data= arrange(DEG_df, p_val_adj, desc(avg_logFC))[1:topFC_labeled,],
# filter(DEG_df, avg_logFC>=1.5)[1:10,],
aes(label=gene), color="black", alpha=.5,
segment.color="black", segment.alpha=.5
) +
# Lines
geom_vline(xintercept= -1.5,lty=4, lwd=.3, alpha=.5) +
geom_vline(xintercept= 1.5,lty=4, lwd=.3, alpha=.5) +
geom_hline(yintercept= -log10(0.05),lty=4, lwd=.3, alpha=.5) +
ggtitle(caption)
print(vol)
}
for (clust in unique(pbmc.markers$cluster)){
cat('\n')
cat("### Cluster ",clust,"\n")
cap <- paste("Cluster",clust,"DEG Table")
DEG_df <- subset(pbmc.markers, cluster==as.character(clust)) %>% arrange(desc(avg_logFC))
volcanoPlot(DEG_df, caption = cap)
createDT(DEG_df, caption = cap)
cat('\n')
}##
## ### Cluster 0
##
##
## ### Cluster 1
##
##
## ### Cluster 2
fp <- FeaturePlot(object = pbmc, features.plot = top1$gene, cols.use = c("grey", "purple"),
reduction.use = "tsne", nCol = nCols, do.return = T)top5 <- pbmc.markers %>% group_by(cluster) %>% top_n(5, avg_logFC)
# setting slim.col.label to TRUE will print just the cluster IDS instead of
# every cell name
DoHeatmap(object = pbmc, genes.use = top5$gene, slim.col.label=T, remove.key=T) %>% ggplotly()RidgePlot(pbmc, features.plot = top1$gene, nCol = nCols, do.sort = F)## Picking joint bandwidth of 0.291
## Picking joint bandwidth of 0.117
## Picking joint bandwidth of 0.0838
markerList <- c("CD14", "FCGR3A")
get_markerDF <- function(pbmc, markerList){
marker.matrix <- pbmc@scale.data[row.names(pbmc@scale.data) %in% markerList, ]
markerMelt <- reshape2:::melt.matrix(marker.matrix)
colnames(markerMelt) <- c("Gene", "Cell", "Expression")
# Fuse metadata
# clusterData <- data.frame(pbmc@ident)
# clusterData$Cell <- row.names(clusterData)
# colnames(clusterData) <- c("Cluster","Cell")
# markerDF <- merge(markerMelt, clusterData, by = "Cell")
metaSelect <- pbmc@meta.data[,c("barcode", "dx", "mut","post_clustering",
"percent.mito","nGene", "nUMI")]
markerDF <- merge(markerMelt,metaSelect, by.x="Cell", by.y="barcode")
return(markerDF)
}
markerDF <- get_markerDF(pbmc, markerList)
createDT(markerDF, caption = "Known Marker Expression")# Explore expression differences between groups
marker_vs_metadata <- function(markerDF, meta_var){
# Create title from ANOVA summary
ANOVAtitle <- function(markerDF, marker){
nTests <- length(unique(markerDF$Gene))
res <- anova(lm(data = subset(markerDF, Gene==marker),
formula = Expression ~ eval(parse(text=meta_var))))
title <-paste(paste("ANOVA (",marker, " vs. ",meta_var, ")", sep=""),
": p=",round(res$`Pr(>F)`,3),
", F=",round(res$`F value`,3),
ifelse(res$`Pr(>F)`<.05/nTests,"(Significant**)",
"(Non-significant)") )
}
title = ""
for (marker in unique(markerDF$Gene) ){
print(marker)
title <- paste(title, "\n", ANOVAtitle(markerDF, marker))
}
ggplot(markerDF, aes(x=eval(parse(text=meta_var)), y=Expression, fill= Gene)) +
geom_boxplot() +
labs(title = title, x=meta_var) +
theme(plot.title = element_text( size=10)) +
scale_fill_manual(values=c("brown", "slategray"))
}marker_vs_metadata(markerDF, "dx")## [1] "CD14"
## [1] "FCGR3A"
marker_vs_metadata(markerDF, "mut") ## [1] "CD14"
## [1] "FCGR3A"
A simplistic way of categorizing cells into CD14++/CD16+ and CD14++/CD16–, is by splitting cells into groups based on whether their expression is higher or lower than the average CD16 expression of all cells.
orig_meta.data <- pbmc@meta.data
################################
avgMarkerExp <-markerDF %>% group_by(Gene) %>% dplyr::summarise(meanExp = mean(Expression))
avgMarkerExp <- setNames(avgMarkerExp$meanExp, avgMarkerExp$Gene)
CD16 <- markerDF[markerDF$Gene=="FCGR3A",]
CD16_group <- ifelse(CD16$Expression >= avgMarkerExp["FCGR3A"],
"CD14++/CD16+", "CD14++/CD16--")
CD16["CD16_group"] <- CD16_group
# Make sure row order is same before putting back into meta.data
metaD <- pbmc@meta.data
newMeta <- merge(metaD, CD16[,c("Cell","CD16_group")], by.x="barcode", by.y="Cell")
row.names(newMeta) <- row.names(metaD)
pbmc <- AddMetaData(pbmc, metadata = newMeta)
# Get proportions of cell types in each cluster
cluster_proportions <- newMeta %>% group_by(CD16_group, post_clustering) %>%
tally() %>%
group_by(CD16_group) %>%
mutate(percentTotal = n/sum(n)*100)
ggplot(cluster_proportions, aes(x=post_clustering, y=percentTotal, fill=CD16_group)) + geom_col(position = "fill") +
labs(title="Proportions of Cell-types per Cluster",
x="Cluster", y="Cell Type / Total Cells") +
scale_fill_manual(values=c("brown", "slategray"))tSNE_metadata_plot("CD16_group")## [1] "t-SNE Metadata plot for CD16_group"
# Show mean exp for each marker
avgMarker <- markerDF %>% group_by(Gene, Cluster) %>% summarise(meanExp = mean(Expression)) ## Error in grouped_df_impl(data, unname(vars), drop): Column `Cluster` is unknown
ggplot(data = avgMarker, aes(x=Gene, y=Cluster, fill=meanExp)) %>% + geom_tile() %>% + scale_fill_distiller(palette="viridis") %>% ggplotly()## Error in ggplot(data = avgMarker, aes(x = Gene, y = Cluster, fill = meanExp)): object 'avgMarker' not found
# Show mean exp for each marker
avgMarker <- markerDF %>% group_by(Gene, dx) %>% summarise(meanExp = mean(Expression))
ggplot(data = avgMarker, aes(x=Gene, y=Cluster, fill=meanExp)) %>% + geom_tile() %>% + scale_fill_distiller(palette="viridis") %>% ggplotly()## Warning in pal_name(palette, type): Unknown palette viridis
## Error in FUN(X[[i]], ...): object 'Cluster' not found
markerMelt <- reshape2::acast(markerDF, Cell~Gene, value.var="Expression", fun.aggregate = mean, drop = F, fill = 0)
#plot_ly( z = markerMelt, y=row.names(markerMelt), z=colnames(markerMelt), type="heatmap")
# dx_colors <- colorRampPalette(brewer.pal(2, "RdBu"))
# mut_colors <- colorRampPalette(brewer.pal(length(unique(pbmc@meta.data$mut)), "Set3"))
Spectral <- grDevices::colorRampPalette(RColorBrewer::brewer.pal(length(unique(pbmc@meta.data$mut)), "Spectral"))
# Spectral <- heatmaply::Spectral(length(unique(pbmc@meta.data$mut)))
heatmaply::heatmaply(markerMelt, key.title="Expression",#plot_method= "plotly",
k_row = length(unique(pbmc.markers)), dendrogram = "row",
showticklabels = c(T, F), xlab = "Known Markers", ylab = "Cells", column_text_angle = 0,
row_side_colors = pbmc@meta.data[,c("dx","mut")], row_side_palette = Spectral
) %>% colorbar(tickfont = list(size = 12), titlefont = list(size = 14), which = 2) %>%
colorbar(tickfont = list(size = 12), titlefont = list(size = 14), which = 1) ## Warning: Didn't find a colorbar to modify.
## Warning: Didn't find a colorbar to modify.
ggplot(data = markerDF, aes(x=post_clustering, y=Expression, fill=Gene)) %>%
+ geom_boxplot(alpha=0.5) %>% + scale_fill_manual(values=c("purple", "turquoise")) # %>% ggplotly() expressionTSNE <- function(pbmc, marker, colors=c("grey", "red")){
FeaturePlot(object = pbmc, features.plot = marker, cols.use = colors,
reduction.use = "tsne", nCol=2, do.return = T, dark.theme = T)[[1]] %>% ggplotly()
}
tp1 <- expressionTSNE(pbmc, markerList[1])tp2 <- expressionTSNE(pbmc, markerList[2], colors=c("grey", "green"))subplot(tp1, tp2)current.cluster.ids <- unique(pbmc.markers$cluster) #c(0, 1, 2, 3, 4, 5, 6, 7)
top1 <- pbmc.markers %>% group_by(cluster) %>% top_n(1, avg_logFC)
new.cluster.ids <- top1$gene #c("CD4 T cells", "CD14+ Monocytes", "B cells", "CD8 T cells", "FCGR3A+ Monocytes", "NK cells", "Dendritic cells", "Megakaryocytes")
pbmc@ident <- plyr::mapvalues(x = pbmc@ident, from = current.cluster.ids, to = new.cluster.ids)
TSNEPlot(object=pbmc, do.label=T, pt.size=0.5, do.return=T) %>% ggplotly()# Available DGE methods:
## "wilcox", "bimod", "roc", "t", "tobit", "poisson", "negbinom", "MAST", "DESeq2"
runDGE <- function(pbmc, meta_var, group1, group2, test.use="wilcox"){
#print(paste("DGE_allCells",meta_var,sep="_"))
pbmc <- SetAllIdent(pbmc, id = meta_var)
pbmc <- StashIdent(pbmc, save.name = meta_var)
DEGs <- FindMarkers(pbmc, ident.1=group1, ident.2=group2, test.use=test.use)
DEGs$gene <- row.names(DEGs)
return(DEGs)
}DEG_df <-runDGE(pbmc, "dx", group1 = "PD", group2="control")
cap = paste("DEGs (All Cells): PD vs. Controls")
createDT(DEG_df, caption = cap)volcanoPlot(DEG_df, caption = cap)DEG_df <-runDGE(pbmc, "mut", "LRRK2", "PD")
cap <- paste("DEGs (All Cells): LRRK2 vs. PD")
createDT(DEG_df, caption = cap)volcanoPlot(DEG_df, caption = cap)DEG_df <-runDGE(pbmc, "CD16_group", "CD14++/CD16+", "CD14++/CD16--")
cap <- paste("DEGs (All Cells): CD14++/CD16+ vs. CD14++/CD16--")
createDT(DEG_df, caption = cap)volcanoPlot(DEG_df, caption = cap)for (clust in unique(pbmc@ident)){
# Subset cells by cluster
pbmc <- SetAllIdent(pbmc, id = "post_clustering")
pbmc_clustSub <- SubsetData(pbmc, ident.use = clust, subset.raw = T)
cap <- paste("Cluster ",clust,": PD vs. Control", sep="")
cat('\n')
cat("### ",cap)
# DGE
DEG_df <-runDGE(pbmc_clustSub, "dx", group1 = "PD", group2="control")
# Show results
volcanoPlot(DEG_df, caption = cap)
createDT(DEG_df, caption = cap)
cat('\n')
} ## Error in WhichCells(object = object, ident = ident.use, ident.remove = ident.remove, : Identity : S100A12 not found.
for (clust in unique(pbmc@ident)){
# Subset cells by cluster
pbmc <- SetAllIdent(pbmc, id = "post_clustering")
pbmc_clustSub <- SubsetData(pbmc, ident.use = clust, subset.raw = T)
cap <- paste("Cluster ",clust,": LRRK2 vs. PD", sep="")
cat('\n')
cat("### ",cap)
# DGE
DEG_df <-runDGE(pbmc_clustSub, "mut", group1 = "LRKK2", group2="PD")
# Show results
volcanoPlot(DEG_df, caption = cap)
createDT(DEG_df, caption = cap)
cat('\n')
} ##
## ### Cluster 1: LRRK2 vs. PD
## Error in WhichCells(object = object, ident = ident.1): Identity : LRKK2 not found.
for (clust in unique(pbmc@ident)){
# Subset cells by cluster
pbmc <- SetAllIdent(pbmc, id = "post_clustering")
pbmc_clustSub <- SubsetData(pbmc, ident.use = clust, subset.raw = T)
cap <- paste("Cluster ",clust,": CD14++/CD16+ vs. CD14++/CD16--", sep="")
cat('\n')
cat("### ",cap)
# DGE
DEG_df <-runDGE(pbmc_clustSub, "CD16_group",
group1 = "CD14++/CD16+", group2="CD14++/CD16--")
# Show results
volcanoPlot(DEG_df, caption = cap)
createDT(DEG_df, caption = cap)
cat('\n')
} ##
## ### Cluster 1: CD14++/CD16+ vs. CD14++/CD16--
##
##
## ### Cluster 0: CD14++/CD16+ vs. CD14++/CD16--
##
##
## ### Cluster 2: CD14++/CD16+ vs. CD14++/CD16--
Further subdivisions within cell types.
If you perturb some of our parameter choices above (for example, setting resolution=0.8 or changing the number of PCs), you might see the CD4 T cells subdivide into two groups. You can explore this subdivision to find markers separating the two T cell subsets. However, before reclustering (which will overwrite object@ident), we can stash our renamed identities to be easily recovered later.
# First lets stash our identities for later
pbmc <- StashIdent(object = pbmc, save.name = "ClusterNames_0.6")
# Note that if you set save.snn=T above, you don't need to recalculate the
# SNN, and can simply put: pbmc <- FindClusters(pbmc,resolution = 0.8)
pbmc <- FindClusters(object = pbmc, reduction.type = "pca", dims.use = 1:10,
resolution = 0.8, print.output = FALSE)## Warning in BuildSNN(object = object, genes.use = genes.use, reduction.type
## = reduction.type, : Build parameters exactly match those of already
## computed and stored SNN. To force recalculation, set force.recalc to TRUE.
## Warning in BuildSNN(object = object, genes.use = genes.use, reduction.type
## = reduction.type, : Build parameters exactly match those of already
## computed and stored SNN. To force recalculation, set force.recalc to TRUE.
# Demonstration of how to plot two tSNE plots side by side, and how to color
# points based on different criteria
plot1 <- TSNEPlot(object = pbmc, do.return = TRUE, no.legend = TRUE, do.label = TRUE, label.size=labSize)
plot2 <- TSNEPlot(object = pbmc, do.return = TRUE, group.by = "ClusterNames_0.6",
no.legend = TRUE, do.label = TRUE, label.size=labSize)
plot_grid(plot1, plot2)# Find discriminating markers
tcell.markers <- FindMarkers(object = pbmc, ident.1 = 0, ident.2 = 1)
# Most of the markers tend to be expressed in C1 (i.e. S100A4). However, we
# can see that CCR7 is upregulated in C0, strongly indicating that we can
# differentiate memory from naive CD4 cells. cols.use demarcates the color
# palette from low to high expression
FeaturePlot(object = pbmc, features.plot = top1$gene, cols.use = c("green", "blue"))pbmc <- SetAllIdent(object = pbmc, id = "ClusterNames_0.6") # Save results for EACH run (in their respective subfolders)
saveRDS(pbmc, file=file.path(params$resultsPath, "cd14-processed.rds") )